function [H,info] = symNMF(A,varargin)        
    par.k = 6;
    par.solver = 'pgd';
    par = process_parameter(par, varargin{:});
    par.m = size(A,1);
    if par.m ~= size(A,2) || ~isequal(A, A')
        error('Need symmetric matrix')
    end
    
    H = rand(par.m,par.k);
    H = eval([par.solver '(A,H,varargin{:})']);
end

function H = pgd(A,H,varargin)
    par.step = 1;
    par.optim_iter = 100;
    par = process_parameter(par, varargin{:});
    
    tAA = mytrace(A,A);
    f = @(H) 0.5*tAA + 0.5*mytrace(H'*H, H'*H)-mytrace(H,A*H);
    step = par.step;
    f0 = f(H);
    iter = 1;
    while step > 1e-12 && iter <= par.optim_iter
        G = 2*H*(H'*H) - 2*A*H;
        H_next = H - step*G;      % gradient
        H_next(H_next < 1e-12) = 0; % project
        f_next = f(H_next);         % calculate f
        if f_next < f0              % line search
            H = H_next;
            f0 = f_next;
            step = par.step;
            iter = iter+1;
            %display(['f0=' num2str(f0)]);
        else
            step = step * 0.5;
        end
    end
end

function t = mytrace(A,B)
    t = sum(sum(A.*B));
end